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Abstract 

Quartic gravity theory is considered with the Einstein-Hilbert La- 



grangean R + aR 2 + bRu U R^ v ', Ruv being Ricci's tensor and R the 



curvature scalar. The parameters a and b are taken of order 1 km 2 . 
Arguments are give which suggest that the effective theory so obtained 
may be a plausible approximation of a viable theory. A numerical in- 
tegration is performed of the field equations for a free neutron gas. 
The result is that the star mass increases with increasing central den- 
sity until about 1 solar mass and then decreases. The baryon number 
increases monotonically, which suggests that the theory allows stars 
in equilibrium with arbitrary baryon number, no matter how large. 
PACS numbers: 04.40.Dg, 04.60.-m, 04.50.Kd 

1 Introduction 

General relativity has passed all observation tests so far, but the real the- 
ory of gravity may well differ significantly from it in strong field regions. In 
fact conceptual difficulties in quantizing Einstein's theory and astrophysical 
observations suggest that general relativity may require modifications. In 
recent years a great effort has been devoted to the study of extended gravity 
theories, mainly with the goal of finding physical explanations to the accel- 
erated expansion of the universe and other astrophysical observations, like 
the flat rotation curves in galaxies pQ, [2]. 

Compact stars are an ideal natural laboratory to look for possible mod- 
ifications of Einsteins theory and their observational signatures. A rather 
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general class of alternative theories of gravity has been considered recently[3] 
to study slowly rotating compact stars with the purpose of investigating con- 
straints on alternative theories. Several studies of compact stars, in particular 
neutron stars, have been made within extended gravity theories[l], [6], [7], 
[5] . There are also theories which prevent the appareance of singularities like 
"gravastars" [S] and Eddington inspired gravity [9]. 

The most popular modification of general relativity, since the early days 
of general relativity [10], derives from an extension of the Einstein-Hilbert 
action of the form 

s = hz I d * x ^( R + F ) + s <**> w 

where k is 8tt times Newton's constant and I use units c = 1 throughout, 
and F is a function of the scalars which may be obtained by combining the 
Riemann tensor, R^xa, and its derivatives, with the metric tensor, g^. In 
particular the theory derived from the choice F(R), where R is the Ricci 
scalar, has been extensively explored under the name of f( R)-gravity [TT] . [12] . 
More general is fourth order gravity [13], which derives from the choice 

F = F (i?, R^R^, R llv \ a R tlu CT ) , (2) 

R^ v being the Ricci tensor. A particular example of eq.flS]) is the quadratic 
Lagrangian which may be written without loss of generality 

F = aR 2 + bR^R^. (3) 

(Riemann square does not appear because it may be eliminated using the 
the Gauss-Bonnett combination 

p2 — p2 A TD D/lf I D T>lLv\a 
K GB = K ~ 4i V^ + ttfiuXa-K , 

which does not contribute to the field equations in a quadratic Lagrangian.) 
The Newtonian limit of the field equations derived from the Lagrangian eq. ([3]) 
has been studied elsewhere [H], [15] . 

In this paper I report on a calculation of neutron stars using the theory 
derived from eq.Q with the particular choice of parameters b = —2a, y/a = 1 
km. That theory is apparently not viable for two reasons. Firstly, in order 
not contradicting Solar System and terrestrial observations the parameters a 
and b should be not greater than a few millimeters. Secondly the weak field 
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limit of the theory should not present ghosts [To]. [T7] . A solution to both 
problems is to assume that eq.flH]) is an approximation, valid for the strong 
fields appearing in neutron stars, of another function F which is extremely 
small in the weak field limit. This would be the case, for instance, if F has 
the form 

F = aR 2 + bR^R^ - c log (l + (a/c)R 2 + (b/c)R^R^) . (4) 

with a ~ 10 6 m 2 ,6 — 2a, c = l/(10 26 m 2 ). Thus eq.()l]) may be approximated 
by 

F~±- c {aR 2 + bR^) 2 <10- 32 R, (5) 

for the Solar System and the relative error due to the terms neglected in 
going from eq.fll]) to eq.([5]) is smaller than 10~ 12 . I have taken into account 
that R 2 ~ R I1V BP V ~ ikp) 2 and that the typical density p ~ 10 4 kg/m 3 . 
The inequality in ([5]) shows that the correction to GR due to the function 
F, eq. 01]) , is neglible in Solar System or terrestrial calculations. Also the 
problem of the ghosts in the weak field limit disappears with that choice. 
Indeed the theory is fine in the context of low-energy effective actions because 
the contribution of R^R^ is so small that it never dominates the dynamics 
of the background. On the other hand the R 2 term does not introduce extra 
graviton modes. 

In contrast in neutron stars, where p ~ 10 18 kg/m 3 , the latter (logaritmic) 
term of eq.fll]) is about 10~ 12 times the former terms and it may be ignored 
in the calculation. 

2 Field equations 

The tensor field equation derived from the functional eqs.flJJ) and the latter 
eq. ([5]) may be taken from the literature [18], [19]. I shall write it in terms of 
the Einstein tensor, G^ u , rather than the Ricci tensor, R^, and in a form 
that looks like the standard Einstein equation of general relativity eq.(]6]). 
That is 

= Rp> - l^uR = k (T™ at + T$) , (6) 
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kT'l = -(2a + b) [V,V U G - g, u UG] - 2 (a + b) 



-GGfiu + -QfiuG 2 



-b 



2G^ l G au — -g^GxcrG a — Vo-V^G^ — V ct V m (j£ + DG^ 



(7) 



We are interested in static problems of spherical symmetry and will use 
the standard metric 



ds 2 = - exp (0 (r)) dt 2 + exp (a (r)) dr 2 + r W. (8) 

Thus G^ u (r) and T™ ai (r) have 3 independent components each, so that 
including a (r) and (3 (r) there are 8 variables. On the other hand there 
are 8 equations, namely 3 eqsflTj) , 3 more equations giving the independent 
components of G^ u in terms of a and j3 and 2 equations of state relating the 3 
independent components of T™ at . I shall assume local isotropy for matter so 
that one of the latter will be the equality T™ at = T^ at (= T^ at in spherical 
symmetry.) In principle the remaining 7 coupled non-linear equations may 
be solved exactly by numerical methods, as will be explained in Section 4. 

Before proceeding, a note about the signs convention is in order. As is well 
known different authors use different signs in the definition of the relevant 
quantities. Here I shall make a choice which essentially agrees with the one 
of Ref . [H] . It may be summarized as follows 

goo = - exp /?, Gf, u = R^p - -g^R = kT^, T ° = -p. (9) 

After that I shall write the three independent components of eq. (Causing the 
notation 

T ° = -p,Tt=p,T 2 = q,T» = T=p + 2q-p, 

{T m at) o Pmati ^~^mat)\ \} mal ) 2 x^rnatj^ Pmat- (10) 

In the following I will name p, p and q the total density, radial pressure and 
transverse pressure respectively, whilst p mat and p mat will be named matter 
density and pressure respectively (remember that we assume local isotropy 
for matter, that is the equality of radial and transverse matter pressures.) 
The differences p — p mat ,p — Pmat and q — p m at will be named effective density, 
radial pressure and transverse pressure respectively. 
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After some algebra I get for the components of the tensor eq.(!7j) 



Pmat 



-p + (2a + b)e 
+&exp(— a) 
+b 



d 2 T 
dr 2 



2 1 A dT 

a — — 

r 2 J dr 

. / Qi an 



+ (a + b)k(^T 2 + 2Tp) 



■— {q-p)+[~oc'P'-P"- — )H> - ■//) 
r \ 2 r 



2(3' 



Ap + 2kp 2 - -k[p 2 +p 2 + 2q 2 } 



(11) 



Pmat 



p -(2a + b)e~ a ( - r + M ^ + (a + b)k(±T 2 - 2Tp) 



Ap + 2kp 2 - -k [p 2 +p 2 + 2q 2 ] 



+b 

+6exp(— a 



2a' 4 



+ -o)(q-p) + [ + P (p + p) Q 2 ) 



Pmat 



q-(2a + b)e 
+ (a + b)k( l -T 2 -2Tq) + b 



d 2 T / 1 1 1 A dT_ 
dr 2 + \ r + 2 1 2 a j dr 



Aq + 2kq 2 - -k [p 2 + p 2 + 2q 2 ] 



+bexp(—a) 



_a' (3' 2 



(3' 

(q-p) + —(p + p) 



Addition of these 3 equations gives the trace equation, that is 

Tmat = Spmat ~ Pmat = T - (6a + 2b) AT, 

where A is the Laplacean operator in curved space-time 



A = exp(— a) 



d 



1 



dr 2 + [r + 2 13 2 a ')dr 



d 



(13) 



(14) 



(15) 



The quantities GJ^ are related to the metric coefficients a and (3 and their 
derivatives (see e.g. [20] ) , hence to p,p and q, that is 

m + Airr 3 p 



exp(— a) 
0" 



^ 2m a' m — 4irpr 3 , 

r ' 2 r 2 — 2mr ' r 2 — 2mr ' 

g 7rr 2 _|_ r p _|_ gy-j ^ ^ m _|_ 4^3^ (y — m 



4nr 3 p) 



2mr 



(r 2 — 2mrY 



(16) 



where I have used units k = 8tt, c = 1 and the radial derivative of a (/?') is 
labelled a' {(3"). The mass parameter m is defined by 

pi- 
rn = I Attx 2 p{x)dx. (17) 
Jo 

The condition that Einstein tensor, G^ u , is divergence free leads to the hy- 
drostatic equilibrium equation, that is 

* = 2h-rt _W (18) 
ar r 2 



3 Application to neutron stars 

For neutron stars, when are quadratic gravity corrections relevant?. In order 
to answer this question we should estimate the conditions where T^, eq.(J7]) , 
is comparable toT™ at - Terms with derivatives are of order 

aDG ~ (cl/Rq)G, 

Rq being the radius of the hypothetical star. Thus these terms are relevant if 
the dimensionless quantity a/R% is of order unity, which implies that a and 
b should be of order the star radius, that is a few kilometers. Terms without 
derivatives are of order 

aG 2 ~ (akp/c 2 )G, 

similar to those with derivatives. 

In order to solve eqs.(|TTT) to ( |T8i) we need an equation of state (eos), that 
is a relation between p mat and p mat , appropriate for a system of neutrons. For 
the calculation here reported I shall choose the eos of a free (non-interacting) 
neutron gas. In order to make easier the rather involved numerical integration 
of the equations, I will simplify the said eos writing 

Pmat = ZPmat + Cp^ t , C = 2.34, (19) 

where p mat and p mat are in units of 7.2 x 10 18 kg m~ 3 . This equation is correct 
in the limit of high density, where p mat ~ 3p mat , and has the same dependence 
Pmat oc Pmat as the eos of the free neutron gas in the nonrelativistic limit of 
low density. The constant C is so chosen that we get the same result as 
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Oppenheimer and Volkoff[21J for the maximum mass stable star in their 
general relativistic calculation. 

A relevant quantity is the baryon number of the star, N, which may be 
calculated from the baryon number density n(r) via 



in our units. A relation between the number density and the matter density 
(or pressure) may be got from the solution of the equation 



which follows from the first law of thermodynamics. Inserting eq.( !T9|) here 
we get a differential equation which may be easily solved with the condition 
Pmat/ n — ^ Ai for p — )• 0, /i being the neutron mass. I get 



where the unit of baryon number density is fi 1 7.2 x 10 18 kg m 3 . 

4 Neutron stars in extended gravity 

In order to derive the structure of neutron stars in generalized f(R) gravity 
theory, as defined by eqs. ([7]) , we should solve the coupled eqs. (TTTj) to (1201 
plus the hydrostatic equilibium eq.( TT8|) with our choice of the parameters a 
and b, namely b = —2a, yfa = 0.96 km. This choose of a and b makes the 
calculation specially simple. 

We need just 3 amongst the 4 eqs. ( Till to (|T4"1) , because only 3 are inde- 
pendent. I choose eqs.f fT4"|) ,the difference eq.f fTB"]) minus eq.f fT2|) , and eq.f fT2|) . 
which may be rewritten 




(20) 



Pmat Tl 



dp-mat 

dn 




(21) 



dT _ dT' 

dr ' dr 



( 




) 



dT 
dr 



T — T,, 



2a 



mat 



(22) 
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— = h',h = q-p, 



dh 

dr 
dh!_ 
dr 



2 1 1 
- + -f3' - -a' ) h! + exp a 
r 2 2 



2a 



+ kTh - 2k(h + 2p)h 



+ — - — ) h + (— - (3" + -a'(3')(p + p) 



(23) 



Pmat 



p + ak(2Tp - -T 2 - 3p 2 + p 2 + 2q 2 ) 



— 2aexp(— a) 



Ap 



2a' 4 
— + ^ 



h+ [f3" --a' f3' ) (p + p) 



(24) 



where the Laplacean operator A was defined in eq. f|T5l) . Finally we need the 
hydrostatic equilibrium eq.f fT8|) . 

The numerical calculation goes as follows. From the values of all vari- 
ables at a given radial coordinate r, integration of the linear differential 
eqs.( )22|) , (123]) and (fT8|) , taking eq.( lT7|) into account, provides the values of 
m, T, T', h, h' and p at r + dr. Hence the relation (see eq. (|T0|) ) 

p = p + 2q - T = 2h + 3p - T, 

gives p (r + dr) , whence eq. f|T8l) gives p' (r + (ir) which allows obtaining p' (r + <ir) 
After that we have all quantities needed to get d 2 p/dr 2 from the derivative 
of eq. ffTSl) , that is 



d 2 p 
dr 2 



2h ' 2h 1 on , \ 1 o, , / / 



Hence we get p mat from eq.( |24]) taking eqs. (|T5|) and (|T6|) into account, which 
allows obtaining p ma< via the eos eq. (|T9|) , whence T ma< = 3p mat — p ma4 follows 
(remember that we assume local isotropy for matter, that is p mat = q ma t-) 
In this way we obtain all the quantities at r + dr and the process may be 
repeated in order to get the quantities at r + 2dr , and so on. This shows 
that our equations form a consistent system. 

As initial conditions for the differential equations we need the values of the 
following variables at the origin: T (0) , h (0) ,p (0) , T' (0) , h' (0) . The latter 
2 should be taken equal to zero in order that there is no singularity, and 



8 



h (0) = because there is no distinction between radial, p, and transverse 
pressure, q at the origin. We are left with just two free parameters, namely 
p (0) and T(0), but there is a constraint, that is the condition that T — > 
for r — ?• oo. Indeed the matter density and pressure are positive within the 
star, so that p mat (r) = for any r > R, R being the star radius (incidentally, 
there is some contribution to the star mass outside the star surface due to 
the effective density.) For r > R(r) the quantity T (r) (and the density p (r)) 
should decrease rapidly with increasing r. As a consequence only the value 
of p (0) may be chosen at will, whilst the value of T(0) should be so chosen as 
to guarantee the rapid decrease of T(r) for r > R. Consequently I have been 
lead to perform the integration several times for each choice of p(0), with 
a different value of T (0) each time, until I got a value of T (r) sufficiently 
small for large enough r (that is greater then the star radius). This procedure 
presents the practical difficulty that requires a fine tuning of T(0) due to the 
fact that for large r > R the solution of eqs. f j22|) is approximately of the form 

™/ x A f r \ B ( r \ 
T (0) ~ — exp —== H exp — . 

r \^2a) r V Via/ 

Thus the parameter A should be very accurately nil in order that the first 
term does not surpasses the second one at large r. This is specially so if 
the parameter a is small, and this is why I have chosen to study the case 
of a relatively large value of a. Also in order to alleviate the problem I have 
substituted a differential equation for a new variable / for the eqs. fl22|) where 

Thus the condition that / remains bounded for r — > oo replaces the stronger 
condition that T — > and the numerical procedure is less unstable. 

In summary we obtain a one-parameter family of equilibrium stars, one 
for each value of the central total pressure p (0). Table 1 reports the results 
of the calculation. As in the standard (GR) theory of neutron stars|21j the 
radius decreases with increasing central density, whilst the mass increases 
until a maximum value and then decreases. Therefore our theory also pre- 
dicts a maximum mass for equilibrium neutron stars. However there is a 
dramatic difference in the behaviour of the baryon number, which here is 
always increasing with increasing central density. Of course in stars with 
very large central density, matter will not be in the form of neutrons but 
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will consists of a mixture of different particles but I will assume that the 
total baryon number is well defined. Although I have not made a rigorous 
proof, the results of the calculation suggest that there may be equilibrium 
configurations of neutron stars for any baryon number no matter how large. 
A consequence of the strong increase of the baryon number with a decrease of 
the mass implies that the binding energy becomes very large, about 90% of 
the mass in the stars with the highest central density amongst those studied 
here. 

Table 1 also shows that both the baryon number density, n, and the 
matter density, p mat , become very large for moderately large central total 
density. This implies that the effective density, p e ^ = p — p mat is negative in 
the central region of the star although becoming positive near and beyond 
the surface. However neither p e ^ nor p mat have a real physical meaning, 
only the total density p being meaningful, and it remains always positive. 
A similar thing happens with the pressure. The surface relative red shift is 
higher in our theory than in the standard (GR) theory, but the difference is 
not dramatic. 

Table 1. Our calculation. Central total pressure, p(0), central total 
density, p (0) , and central matter density, p mat (0), are in units p c = 7.2 x 10 18 
kg/m 3 . Central baryon number density, n(0), in units p c /p, p being the 
neutron mass. Star radius, R, is in kilometers, mass, M, in solar masses 
and baryon number, N, in units of solar baryon number. I report also the 
dimensionless fractional surface red shift, A A/A = 1/ a/1 — 2M/R — 1, and 
percent binding energy, BE = 100(iV — M)/N . An expressions like 1.6E2 
means 1.6 x 10 2 



p(0) 


0.01 


0.1 


1 


10 


100 


1000 


10000 


P(0) 


0.18 


0.82 


4.5 


34 


3.LE72 


3.0E3 


3.0E5 


Pmat (0) 


1.6E2 


2.5£3 


4.3E4 


7.8E5 


1.QE7 


3.2E8 


6.3E9 


n(0) 


56 


A.5E2 


3.7E3 


3.3EA 


3.2Eh 


3.0E6 


2.8£7 


R 


10.7 


6.7 


4.0 


2.7 


2.1 


2.2 


2.2 


M 


0.67 


0.80 


0.60 


0.39 


0.264 


0.268 


0.292 


N 


0.73 


0.94 


1.03 


1.13 


1.44 


2.00 


2.63 


BE 


8.9% 


15% 


41% 


65% 


82% 


87% 


89% 


AA/A 


0.106 


0.22 


0.34 


0.31 


0.23 


0.23 


0.26 
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5 Discussion 



The calculations of this paper show that, if there are corrections to general 
relativiy of the form of eqs.fJJJ) andQ, then the structure of neutron stars 
would be dramatically different from the one predicted by the general rel- 
ativity. In particular a new scenario would emerge for the evolution of the 
central region of massive white dwarfs stars after the supernova explosion. 
Indeed the said central region might contract strongly by emitting an amount 
of energy corresponding to a very large fraction of the total mass. The fi- 
nal result will be a neutron star with a mass maybe surpassing the believed 
(Oppenheimer-Volkoff) limit. It is not possible to know how large is the 
new mass limit until calculations with more realistic equations of state are 
performed. In addition the predictions of the theory here considered may be 
quite different for other choices of the parameters a and b. 

6 Appendix. Neutron stars in general rela- 
tivity 

For the sake of comparison with the results of our calculation using eqs.f fTT]) 
to (|2ip , I summarize in Table 2 the results of a calculation similar to the 
one performed by the Oppenheimer and Volkoff calculation[21J using general 
relativity. It corresponds to taking a = b = in eqs.f fTT]) to ( !T4|) . so that 
P = PmatiV = Pmat, and I use the eos eq. (|T9|) . I have extended the calculation 
to quite high values of the central pressure because for those values the 
corrections to GR in our generalized f{R) gravity theory are most relevant. 

Table 2. General relativistic calculation. Central pressure, p(0), and 
central density, p (0) , are in units 7.2 x 10 18 kg to -3 , star radius, R, in 
kilometers, mass, M, in solar masses and baryon number, N, in units of 
solar baryon number. I report also the percent binding energy, BE, defined 

by the ratio 100(iV — M)/N and the fractional surface red shift, AA/A = 

l/y/l-2M/R-l . 
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p(0) 


0.01 


0.04 


0.1 


1 


10 


100 


1000 


P(0) 


0.18 


0.46 


0.89 


5.3 


39 


3.4 x 10 2 


3.1 x 10 3 


i? 


11.9 


9.6 


8.3 


5.8 


5.2 


6.6 


6.6 


M 


0.67 


0.72 


0.70 


0.55 


0.39 


0.40 


0.44 


N 


0.69 


0.74 


0.73 


0.55 


0.36 


0.37 


0.42 


BE 


2.9% 


3.4% 


3.4% 


-0.73% 


-8.1% 


-8.0% 


-5.7% 


AA/A 


0.094 


0.13 


0.15 


0.18 


0.13 


0.11 


0.12 



Table 2 shows that both the mass, M, and the baryon number, N, in- 
crease with increasing central density until p (0) ~ 0.46, where M ~ 0.72 
M Q (the OV mass limit), and both M and N decrease after that. There are 
no equilibrium configurations, either stable or unstable, with baryon number 
above N ~ 0.74. Actually for every baryon number N < 0.74 there are two 
equilibrium configurations, one of them with p (0) < 0.46 and another one 
with p (0) > 0.46, the latter having higher mass than the former. Further- 
more, as is shown in Table 2, stars with large central density have a negative 
binding energy and therefore cannot be stable. 
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